https://femiguez.github.io/nlraa-docs/nlraa-Oddi-LFMC.html

https://onlinelibrary.wiley.com/doi/pdfdirect/10.1002/ece3.5543

#install.packages("emmeans")

library(nlme)
library(nlraa)
library(ggplot2)
library(emmeans)   # emmeans::contrast(), emmeans::emmeans()

Steps in fitting a nonlinear mixed model

Part I.

data(lfmc)
lfmc.gd <- groupedData(lfmc ~ time | group, data = lfmc, order.groups = FALSE) 
lfmc.gd$plot = with(lfmc.gd, factor(plot, levels =c("4", "5", "6", "1", "2", "3")))
head(lfmc.gd)
Grouped Data: lfmc ~ time | group
         leaf.type time plot site  lfmc     group
1      M. spinosum    1    1    P 238.4 SM plot 1
2      M. spinosum    1    1    P 269.8 SM plot 1
3      M. spinosum    1    1    P 325.2 SM plot 1
4 S. bracteolactus    1    1    P 326.7 SS plot 1
5 S. bracteolactus    1    1    P 257.3 SS plot 1
6 S. bracteolactus    1    1    P 279.3 SS plot 1
ggplot(lfmc, aes(time, lfmc)) +
  geom_point() +
  geom_smooth(method = "loess") +
  facet_wrap(~leaf.type) +
  xlab("Time (days)") +
  ylab("LMFC (%)")
`geom_smooth()` using formula 'y ~ x'

NONLINEAR MIXED-EFFECTS MODEL (M1)

nlsL <- nlsList(lfmc ~ SSdlf(time, A, w, m, s), data = lfmc.gd)
Warning: 3 errors caught in nls(model, data = data, control = controlvals).  The error messages and their frequencies are

step factor 0.000488281 reduced below 'minFactor' of 0.000976562 
                                                               1 
                     number of iterations exceeded maximum of 50 
                                                               2 
coef(nlsL)
plot(intervals(nlsL))

Random-Effects Mode

nl.re <- nlme(nlsL, control = nlmeControl(maxIter = 5000, msMaxIter = 1500))
Warning in (function (model, data = sys.frame(sys.parent()), fixed, random,  :
  Iteration 1, LME step: nlminb() did not converge (code = 1). PORT message: function evaluation limit reached without convergence (9)
nl.re.1 <- update(nl.re, random = pdDiag(A + w + m + s ~ 1))
fxf <- fixef(nl.re.1)
fxf
        A         w         m         s 
193.45254  21.45745  31.48273 -22.59708 
# nl.me <- update(nl.re.1, fixed = list(A + w + m + s ~ leaf.type),
#                 start = c(A = fxf[1],0,0,0,
#                           w = fxf[2],0,0,0,
#                           m = fxf[3],0,0,0,
#                           s = fxf[4],0,0,0))
nl.re.2 <- update(nl.re.1, random = pdDiag(A + w + m ~ 1))
fxf <- fixef(nl.re.2)
fxf
        A         w         m         s 
193.46654  21.43199  31.48852 -22.60602 
nl.me.1 <- update(
  nl.re.2, 
  fixed = list(A + w + m + s ~ leaf.type),  
  start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3],0,0,0, s=fxf[4],0,0,0))
nl.me.1
Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1034.375
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E 
                  23.135414                   53.294908 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                 296.620596                  263.313118 
              w.(Intercept)          w.leaf.typeGrass E 
                  50.291577                  -34.313924 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                 -26.732287                    2.135665 
              m.(Intercept)          m.leaf.typeGrass E 
                  51.533983                  -21.920404 
     m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                 -20.647403                  -18.170359 
              s.(Intercept)          s.leaf.typeGrass E 
                  17.605794                  -25.926031 
     s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                 -43.529450                  -33.666226 

Random effects:
 Formula: list(A ~ 1, w ~ 1, m ~ 1)
 Level: group
 Structure: Diagonal
        A.(Intercept) w.(Intercept) m.(Intercept) Residual
StdDev:      8.112813   0.001245815      2.468334 15.38607

Number of Observations: 247
Number of Groups: 12 
nl.re.2
Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1070.648
  Fixed: list(A ~ 1, w ~ 1, m ~ 1, s ~ 1) 
        A         w         m         s 
193.46654  21.43199  31.48852 -22.60602 

Random effects:
 Formula: list(A ~ 1, w ~ 1, m ~ 1)
 Level: group
 Structure: Diagonal
               A        w        m Residual
StdDev: 119.3792 10.94165 8.927279 15.26574

Number of Observations: 247
Number of Groups: 12 
AIC(nl.re.2, nl.me.1)
IC_tab(nl.re.2, nl.me.1)
nl.me.2 <- update(nl.me.1, random = pdDiag(A + m ~ 1))
IC_tab(nl.me.1, nl.me.2)
nl.me.3 <- update(nl.me.2, random = pdDiag(A ~ 1)) 
IC_tab(nl.me.2, nl.me.3)
plot(nl.me.3)

rse <- residuals(nl.me.3) / sigma(nl.me.3)
hist(residuals(nl.me.3, type = "normalized"))

#par(pty = "s")
qqnorm(rse)
qqline(rse)

nl.me.3.vm <- update(nl.me.3, weights = varIdent(form = ~1 | leaf.type))
IC_tab(nl.me.3, nl.me.3.vm)
plot(nl.me.3.vm)

hist(residuals(nl.me.3.vm, type = "normalized"), main = "", xlab = "Residuals")

qqnorm(resid(nl.me.3.vm, type = "normalized"))
qqline(resid(nl.me.3.vm, type = "normalized"))

plot(ACF(nl.me.3.vm), alpha = 0.01)

Evaluation of Fixed Effects

nl.me.ml <- update(nl.me.3.vm, method = "ML")
nl.me.ml
Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -965.5628
  Fixed: list(A + w + m + s ~ leaf.type) 
              A.(Intercept)          A.leaf.typeGrass E 
                 50.1522850                  28.9553860 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                346.5816365                 231.1569668 
              w.(Intercept)          w.leaf.typeGrass E 
                 24.6693608                  -9.4837876 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                -35.8327959                  29.7533014 
              m.(Intercept)          m.leaf.typeGrass E 
                 48.9613646                 -19.9780200 
     m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                -26.2920797                 -14.9909275 
              s.(Intercept)          s.leaf.typeGrass E 
                -16.2294126                   6.3520529 
     s.leaf.typeM. spinosum s.leaf.typeS. bracteolactus 
                -21.0664461                   0.9389322 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.412262 6.824811

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8547836        3.1456077        3.0262123 
Number of Observations: 247
Number of Groups: 12 
nl.me.ml.1 <- update(nl.me.ml,
  fixed = list(A + w + m ~ leaf.type, s ~ 1),
  start = c(A = fxf[1], 0, 0, 0, w = fxf[2], 0, 0, 0, m = fxf[3], 0, 0, 0, s = fxf[4])
)
nl.me.ml.1
Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -969.2617
  Fixed: list(A + w + m ~ leaf.type, s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E 
                   49.94562                    39.60979 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  219.45652                   232.11909 
              w.(Intercept)          w.leaf.typeGrass E 
                   25.24501                   -14.34293 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   35.36121                    28.70388 
              m.(Intercept)          m.leaf.typeGrass E 
                   48.30390                   -21.59477 
     m.leaf.typeM. spinosum m.leaf.typeS. bracteolactus 
                  -15.57765                   -14.37454 
                          s 
                  -15.44291 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.612472 6.825604

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8714808        3.2750695        3.0238813 
Number of Observations: 247
Number of Groups: 12 
IC_tab(nl.me.ml, nl.me.ml.1) 
nl.me.ml.2 <- update(nl.me.ml, fixed = list(A + w ~ leaf.type, m + s ~ 1),  
                     start = c(A=fxf[1],0,0,0, w=fxf[2],0,0,0, m=fxf[3], s=fxf[4])) 
nl.me.ml.2
Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -972.7465
  Fixed: list(A + w ~ leaf.type, m + s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E 
                   54.25730                    30.93396 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  223.27914                   240.31543 
              w.(Intercept)          w.leaf.typeGrass E 
                   29.05803                   -20.68663 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   31.68646                    26.99003 
                          m                           s 
                   30.91655                   -16.15435 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.216113 7.028345

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8820413        3.1693967        2.9643305 
Number of Observations: 247
Number of Groups: 12 
nl.me.ml.3 <- update(nl.me.ml, fixed = list(A ~ leaf.type, w + m + s  ~ 1),  
                     start = c(A=fxf[1],0,0,0, w=fxf[2], m=fxf[3], s=fxf[4])) 
nl.me.ml.3
Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -1017.91
  Fixed: list(A ~ leaf.type, w + m + s ~ 1) 
              A.(Intercept)          A.leaf.typeGrass E 
                   63.92214                    14.27256 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  251.99464                   265.42794 
                          w                           m 
                   17.33007                    32.65667 
                          s 
                  -25.16304 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:       8.13893 7.795662

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
        1.000000         1.584893         2.813972         2.731419 
Number of Observations: 247
Number of Groups: 12 
nl.me.ml.4 <- update(nl.me.ml, fixed = list(w ~ leaf.type, A + m + s  ~ 1),  
                     start = c(A=fxf[1], w=fxf[2],0,0,0, m=fxf[3], s=fxf[4])) 
nl.me.ml.4
Nonlinear mixed-effects model fit by maximum likelihood
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 
  Log-likelihood: -999.8224
  Fixed: list(w ~ leaf.type, A + m + s ~ 1) 
              w.(Intercept)          w.leaf.typeGrass E 
                   29.16028                   -20.38030 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   33.19422                    28.72473 
                          A                           m 
                  175.78134                    31.08599 
                          s 
                  -15.70300 

Random effects:
 Formula: A ~ 1 | group
               A Residual
StdDev: 108.0655 7.042353

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8770644        3.1888994        2.9304282 
Number of Observations: 247
Number of Groups: 12 
IC_tab(nl.me.ml.1, nl.me.ml.2, nl.me.ml.3, nl.me.ml.4)
M1 <- update(nl.me.ml.2, method = "REML")
summary(M1)
Nonlinear mixed-effects model fit by REML
  Model: lfmc ~ SSdlf(time, A, w, m, s) 
  Data: lfmc.gd 

Random effects:
 Formula: A ~ 1 | group
        A.(Intercept) Residual
StdDev:      9.408521 7.162899

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W          Grass E      M. spinosum S. bracteolactus 
       1.0000000        0.8851762        3.1796110        2.9647390 
Fixed effects:  list(A + w ~ leaf.type, m + s ~ 1) 
 Correlation: 
                            A.(In) A.l.GE A..M.s A..S.b w.(In) w.l.GE w..M.s
A.leaf.typeGrass E          -0.527                                          
A.leaf.typeM. spinosum      -0.146  0.535                                   
A.leaf.typeS. bracteolactus -0.115  0.537  0.746                            
w.(Intercept)               -0.217 -0.034 -0.200 -0.216                     
w.leaf.typeGrass E          -0.035 -0.283 -0.382 -0.399 -0.258              
w.leaf.typeM. spinosum      -0.118 -0.227 -0.547 -0.454  0.157  0.573       
w.leaf.typeS. bracteolactus -0.129 -0.241 -0.462 -0.575  0.188  0.601  0.640
m                           -0.217 -0.324 -0.633 -0.666  0.092  0.145  0.161
s                           -0.246 -0.354 -0.710 -0.748  0.416  0.575  0.702
                            w..S.b m     
A.leaf.typeGrass E                       
A.leaf.typeM. spinosum                   
A.leaf.typeS. bracteolactus              
w.(Intercept)                            
w.leaf.typeGrass E                       
w.leaf.typeM. spinosum                   
w.leaf.typeS. bracteolactus              
m                            0.172       
s                            0.752  0.544

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.27999716 -0.59657157 -0.02042043  0.57273326  3.12017921 

Number of Observations: 247
Number of Groups: 12 
fixef(M1)
              A.(Intercept)          A.leaf.typeGrass E 
                   54.25350                    30.92293 
     A.leaf.typeM. spinosum A.leaf.typeS. bracteolactus 
                  223.24504                   240.27654 
              w.(Intercept)          w.leaf.typeGrass E 
                   29.05581                   -20.68938 
     w.leaf.typeM. spinosum w.leaf.typeS. bracteolactus 
                   31.67297                    26.97508 
                          m                           s 
                   30.92881                   -16.15406 
ranef(M1)
          A.(Intercept)
GW plot 4     -1.915443
GW plot 5     -1.923198
GW plot 6      3.838641
GE plot 1     15.652459
SM plot 1      6.967303
SS plot 1      9.166925
GE plot 2    -11.066197
SM plot 2     -5.358910
SS plot 2      2.442781
GE plot 3     -4.586263
SM plot 3     -1.608393
SS plot 3    -11.609707
intervals(M1)
Approximate 95% confidence intervals

 Fixed effects:
                                lower      est.     upper
A.(Intercept)                42.45370  54.25350  66.05331
A.leaf.typeGrass E           13.51301  30.92293  48.33286
A.leaf.typeM. spinosum      191.88318 223.24504 254.60691
A.leaf.typeS. bracteolactus 207.05884 240.27654 273.49423
w.(Intercept)                25.68039  29.05581  32.43122
w.leaf.typeGrass E          -25.79788 -20.68938 -15.58088
w.leaf.typeM. spinosum       16.39521  31.67297  46.95073
w.leaf.typeS. bracteolactus  11.07083  26.97508  42.87934
m                            26.85931  30.92881  34.99831
s                           -20.81510 -16.15406 -11.49302

 Random Effects:
  Level: group 
                     lower     est.    upper
sd(A.(Intercept)) 5.035417 9.408521 17.57953

 Variance function:
                     lower      est.    upper
Grass E          0.6772617 0.8851762 1.156919
M. spinosum      2.4528838 3.1796110 4.121649
S. bracteolactus 2.2864085 2.9647390 3.844316

 Within-group standard error:
   lower     est.    upper 
5.970704 7.162899 8.593144 

Part II

# 2.2.1 - Overall predictions (Table 3) ---- 
Agw <- fixef(M1)[1] # "A" for grasses in the W site (GW) 
Age <- fixef(M1)[1]+fixef(M1)[2] # "A" for grasses in the E site (GE)
Asm <- fixef(M1)[1]+fixef(M1)[3] # "A" for Mullinum spinosum (SM)
Ass <- fixef(M1)[1]+fixef(M1)[4] # "A" for Senecio filaginoides (SS)
Wgw <- fixef(M1)[5] # "w" for grasses in the W site (GW) 
Wge <- fixef(M1)[5]+fixef(M1)[6] # "w" for grasses in the E site (GE) 
Wsm <- fixef(M1)[5]+fixef(M1)[7] # "w" for Mullinum spinosum (SM)
Wss <- fixef(M1)[5]+fixef(M1)[8] # "w" for Senecio filaginoides (SS)
M <- fixef(M1)[9] # "m" for all the leaf types
S <- fixef(M1)[10] # "s" for all the leaf types

GW <- subset(lfmc, leaf.type == "Grass W")
GE <- subset(lfmc, leaf.type == "Grass E")
SM <- subset(lfmc, leaf.type == "M. spinosum")
SS <- subset(lfmc, leaf.type == "S. bracteolactus")

par(mfrow=c(2,2), cex=0.75)
# Grasses W
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=2, lty=1, add=T, col="darkgreen") 
# Grasses E
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age-Wge)/(1 + exp((M-x)/S))+Wge), lwd=2, lty=1, add=T, col="green")
# M. Spinosum
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=2, lty=1, add=T, col="darkorange") 
# S. filaginoides 
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass-Wss)/(1 + exp((M-x)/S))+Wss), lwd=2, lty=1, add=T, col="darkred") 

# 2.2.2 - Predictions at plot level (Table 3) ----
Agw.p4 <- fixef(M1)[1]+ranef(M1)[1,1] # "A" for grasses in the plot 4 (GW) 
Agw.p5 <- fixef(M1)[1]+ranef(M1)[2,1] # "A" for grasses in the plot 5 (GW) 
Agw.p6 <- fixef(M1)[1]+ranef(M1)[3,1] # "A" for grasses in the plot 6 (GW) 
Age.p1 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[4,1] # "A" for grasses in the plot 1 (GW) 
Age.p2 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[7,1] # "A" for grasses in the plot 2 (GW) 
Age.p3 <- fixef(M1)[1]+fixef(M1)[2]+ranef(M1)[10,1] # "A" for grasses in the plot 3 (GW) 
Asm.p1 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[5,1] # "A" for M. spinosum in the plot 1 (SM) 
Asm.p2 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[8,1] # "A" for M. spinosum in the plot 2 (SM)
Asm.p3 <- fixef(M1)[1]+fixef(M1)[3]+ranef(M1)[11,1] # "A" for M. spinosum in the plot 3 (SM)
Ass.p1 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[6,1] # "A" for S. filaginoides in the plot 1 (SM)   
Ass.p2 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[9,1] # "A" for S. filaginoides in the plot 2 (SM) 
Ass.p3 <- fixef(M1)[1]+fixef(M1)[4]+ranef(M1)[12,1] # "A" for S. filaginoides in the plot 3 (SM)
  
par(mfrow=c(2,2), cex=0.75) # (Figure 4)
# Grasses W site 
plot(GW$time, GW$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses W site", col="darkgreen")
curve(((Agw.p4-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 4
curve(((Agw.p5-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 5
curve(((Agw.p6-Wgw)/(1 + exp((M-x)/S))+Wgw), lwd=1, lty=2, add=T, col="darkgreen") # plot 6
# Grasses E site
plot(GE$time, GE$lfmc, xlab="", ylab="LFMC (%)", ylim=c(0,100), main="Grasses E site", col="green")
curve(((Age.p1-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 1
curve(((Age.p2-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 2
curve(((Age.p3-Wge)/(1 + exp((M-x)/S))+Wge), lwd=1, lty=2, add=T, col="green") # plot 3
# M. Spinosum (E site) 
plot(SM$time, SM$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="M. spinosum", font.main=4, col="darkorange")
curve(((Asm.p1-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 1
curve(((Asm.p2-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 2
curve(((Asm.p3-Wsm)/(1 + exp((M-x)/S))+Wsm), lwd=1, lty=2, add=T, col="darkorange") # plot 3
# S. filaginoides (E site)
plot(SS$time, SS$lfmc, xlab="Time (days)", ylab="LFMC (%)", ylim=c(0,350), main="S. bracteolactus", font.main=4, col="darkred")
curve(((Ass.p1-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 1
curve(((Ass.p2-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 2
curve(((Ass.p3-Wss)/(1 + exp((M-x)/S))+Wss), lwd=1, lty=2, add=T, col="darkred") # plot 3

# 2.2.3 - Derivatives (drying speed) ---- 
lfmc.GW <- expression((Agw-Wgw)/(1 + exp((M-x)/S))+Wgw) # LFMC dynamics for grasses in the W site
ds.GW <- deriv(lfmc.GW, "x", function.arg = TRUE) # Drying speed for grasses in the W site

lfmc.GE <- expression((Age-Wge)/(1 + exp((M-x)/S))+Wge) # LFMC dynamics for grasses in the E site
ds.GE <- deriv(lfmc.GE, "x", function.arg = TRUE) # Drying speed for grasses in the E site

lfmc.SM <- expression((Asm-Wsm)/(1 + exp((M-x)/S))+Wsm) # LFMC dynamics for M. spinosum
ds.SM <- deriv(lfmc.SM, "x", function.arg = TRUE) # Drying speed for M. spinosum

lfmc.SS <- expression((Ass-Wss)/(1 + exp((M-x)/S))+Wss) # LFMC dynamics for S. filaginoides
ds.SS <- deriv(lfmc.SS, "x", function.arg = TRUE) # # Drying speed for S. filaginoides

xvec <- seq(0, 100, length = 1000)
y.ds.GW <- ds.GW(xvec)
y.ds.GE <- ds.GE(xvec)
y.ds.SM <- ds.SM(xvec)
y.ds.SS <- ds.SS(xvec)

par(mfrow=c(2,2), cex=0.75) # (Figure 4)
plot(xvec, y.ds.GW, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the W site")
lines(xvec, -1*(attr(y.ds.GW, "grad")), lwd=2, lty=1, col="darkgreen")
plot(xvec, y.ds.GE, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="Grasses in the E site")
lines(xvec, -1*(attr(y.ds.GE, "grad")), lwd=2, lty=1, col="green")
plot(xvec, y.ds.SM, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="M. spinosum", font.main=4)
lines(xvec, -1*(attr(y.ds.SM, "grad")), lwd=2, lty=1, col="darkorange")
plot(xvec, y.ds.SS, xlim=c(0,80), ylim=c(0,4), 
     ylab="Drying speed", xlab="Time (days)", type="n", main="S. filaginoides", font.main=4)
lines(xvec, -1*(attr(y.ds.SS, "grad")), lwd=2, lty=1, col="darkred")

# 2.3 ---- Testing differences among leaf types ------------------------------------------- 

# 2.3.1 - Parameter "A" ----
contrast(emmeans(M1, ~leaf.type, param = "A"), "pairwise")
 contrast                       estimate    SE  df t.ratio p.value
 Grass W - Grass E                 -30.9  8.84 226  -3.500  0.0031
 Grass W - M. spinosum            -223.2 15.92 226 -14.027  <.0001
 Grass W - S. bracteolactus       -240.3 16.86 226 -14.254  <.0001
 Grass E - M. spinosum            -192.3 13.45 226 -14.294  <.0001
 Grass E - S. bracteolactus       -209.4 14.22 226 -14.718  <.0001
 M. spinosum - S. bracteolactus    -17.0 11.71 226  -1.454  0.4671

P value adjustment: tukey method for comparing a family of 4 estimates 
# 2.3.2 - Parameter "w" ---- 
contrast(emmeans(M1, ~leaf.type, param = "w"), "pairwise")
 contrast                       estimate   SE  df t.ratio p.value
 Grass W - Grass E                  20.7 2.59 226   7.981  <.0001
 Grass W - M. spinosum             -31.7 7.75 226  -4.085  0.0004
 Grass W - S. bracteolactus        -27.0 8.07 226  -3.342  0.0053
 Grass E - M. spinosum             -52.4 6.62 226  -7.914  <.0001
 Grass E - S. bracteolactus        -47.7 6.84 226  -6.974  <.0001
 M. spinosum - S. bracteolactus      4.7 6.72 226   0.699  0.8974

P value adjustment: tukey method for comparing a family of 4 estimates 

# 3.1.1 - Base nonlinear model ----
nl.fe.0 <- nls(lfmc ~ ((A - w) / (1 + exp((m - time)/s))) + w, 
               start = c(A=15, m=30, s=-17, w=5),
               data = lfmc) # here, time is the only predictor variable 

b <- coef(nl.fe.0) # coefficients of the base model

# 3.1.2 - Leaf type fixed effect ----
nl.fe.1 <- nls(lfmc ~ ((A[leaf.type] - w[leaf.type])/(1 + exp((m - time)/s))) + w[leaf.type], 
              start = list(A=rep(b[1], 4), m=25, s=-10, w=rep(b[4], 4)),
              data = lfmc) # and fit is made using the base model's coefficients as start values

b1 <- coef(nl.fe.1) # coefficients of the model with leaf type and time as predictor variables

# 3.1.3 - Plot fixed effect ---- 
nl.fe.2 <- nls(lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group], 
               start = list(A=rep(c(b1[1],b1[2],b1[3],b1[4]),3), m=25, s=-10, w=rep(c(b1[7],b1[8],b1[9],b1[10]),3)),
               data = lfmc) # the model is fit using "b1" as the start values

IC_tab(nl.fe.0, nl.fe.1, nl.fe.2)
# Residual checking:
par(mfrow=c(1,1))
plot(nl.fe.2) # heterocedasticity in the residuals

hist(resid(nl.fe.2), main="", xlab="Residuals") # slightly skew distribution   

qqnorm(resid(nl.fe.2))
qqline(resid(nl.fe.2))

M2 <- nl.fe.2
summary(M2)

Formula: lfmc ~ ((A[group] - w[group])/(1 + exp((m - time)/s))) + w[group]

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
A1    53.683      8.602   6.240 2.20e-09 ***
A2    54.282      8.636   6.285 1.72e-09 ***
A3    61.906      8.881   6.971 3.59e-11 ***
A4   111.561     13.291   8.394 5.65e-15 ***
A5   314.131     24.840  12.646  < 2e-16 ***
A6   333.831     26.640  12.531  < 2e-16 ***
A7    77.465     10.997   7.044 2.34e-11 ***
A8   298.161     24.917  11.966  < 2e-16 ***
A9   321.433     25.765  12.475  < 2e-16 ***
A10   87.553     11.744   7.455 2.01e-12 ***
A11  279.437     20.829  13.416  < 2e-16 ***
A12  292.443     24.178  12.095  < 2e-16 ***
m     30.413      2.896  10.504  < 2e-16 ***
s    -20.695      3.606  -5.739 3.11e-08 ***
w1    28.369      6.537   4.340 2.17e-05 ***
w2    27.480      6.548   4.197 3.93e-05 ***
w3    25.962      6.633   3.914 0.000121 ***
w4     0.890      8.173   0.109 0.913388    
w5    43.524     13.567   3.208 0.001535 ** 
w6    41.223     14.429   2.857 0.004686 ** 
w7     6.091      7.262   0.839 0.402466    
w8    26.610     13.604   1.956 0.051725 .  
w9    39.495     14.009   2.819 0.005252 ** 
w10    2.265      7.551   0.300 0.764455    
w11   66.661     11.478   5.808 2.19e-08 ***
w12   38.222     13.031   2.933 0.003708 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.62 on 221 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 7.269e-06
# 4.1.1 - Base linear mixed-effect model ----
l.me <- lme(lfmc ~ time * leaf.type, random = ~1 | plot, data = lfmc)

# Residual checking:
plot(l.me) # heterocedasticity in the residuals

plot(ACF(l.me, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # temporal correlation is observed at lag 1

# 4.1.2 - Variance modelling ---- 
l.me.vm <- update(l.me, weights = varIdent(form=~ 1 | leaf.type))
    
# Residual checking:
plot(l.me.vm) # variances there seem to be well modeled

AIC(l.me, l.me.vm) # modeling the variance improves the model fit 
# 4.1.3 - Temporal correlation modelling ----  
l.me.vm.tm <- update(l.me.vm, correlation=corARMA(p=1,q=1, form=~1))
plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) # the temporal dependence is modeled

AIC(l.me.vm, l.me.vm.tm) # modeling the temporal correlation improves the model fit
plot(l.me.vm.tm) 

plot(ACF(l.me.vm.tm, resType = "normalized", na.action=na.omit), alpha=0.05, grid=TRUE) 

hist(resid(l.me.vm.tm, type="normalized"), main="", xlab="Residuals") 

qqnorm(resid(l.me.vm.tm, type="normalized"))
qqline(resid(l.me.vm.tm, type="normalized"))

# 4.1.4 - Evaluation of the fixed-effects ---- 
l.me.ml <- update(l.me.vm.tm, method ="ML") # the model is fitted by Maximum likelihood 
l.me.ml.1 <- update(l.me.ml, ~.-time:leaf.type) # the model is reduced removing the interaction "time x leaf type" 
IC_tab(l.me.ml, l.me.ml.1) # the AIC comparison suggests the interaction term to be important
M3 <- l.me.vm.tm 
summary(M3)
Linear mixed-effects model fit by REML
  Data: lfmc 

Random effects:
 Formula: ~1 | plot
        (Intercept) Residual
StdDev:    3.750488 7.842661

Correlation Structure: ARMA(1,1)
 Formula: ~1 | plot 
 Parameter estimate(s):
      Phi1     Theta1 
 0.9076001 -0.7567060 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | leaf.type 
 Parameter estimates:
         Grass W      M. spinosum S. bracteolactus          Grass E 
        1.000000         2.963571         3.123740         1.070897 
Fixed effects:  lfmc ~ time * leaf.type 
 Correlation: 
                               (Intr) time   lf.tGE lf.M.s lf.S.b tm:.GE
time                           -0.557                                   
leaf.typeGrass E               -0.709  0.395                            
leaf.typeM. spinosum           -0.407  0.227  0.654                     
leaf.typeS. bracteolactus      -0.394  0.219  0.653  0.666              
time:leaf.typeGrass E           0.374 -0.671 -0.592 -0.430 -0.429       
time:leaf.typeM. spinosum       0.176 -0.316 -0.299 -0.742 -0.409  0.546
time:leaf.typeS. bracteolactus  0.166 -0.299 -0.318 -0.441 -0.743  0.562
                               t:.M.s
time                                 
leaf.typeGrass E                     
leaf.typeM. spinosum                 
leaf.typeS. bracteolactus            
time:leaf.typeGrass E                
time:leaf.typeM. spinosum            
time:leaf.typeS. bracteolactus  0.567

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.3449313 -0.6438200 -0.1011963  0.6158869  3.3469739 

Number of Observations: 247
Number of Groups: 6 
M4 <- lm(lfmc ~ time * group, data = lfmc)
summary(M4)

Call:
lm(formula = lfmc ~ time * group, data = lfmc)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.673  -7.042   0.466   7.002  66.332 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          48.79275    6.44309   7.573 9.60e-13 ***
time                 -0.24478    0.13305  -1.840  0.06714 .  
groupGW plot 5        0.49817    9.11190   0.055  0.95645    
groupGW plot 6        6.14984    9.11190   0.675  0.50042    
groupGE plot 1       40.32675    9.49655   4.246 3.19e-05 ***
groupSM plot 1      212.68278    9.11190  23.341  < 2e-16 ***
groupSS plot 1      227.42468    9.11190  24.959  < 2e-16 ***
groupGE plot 2       14.39850    9.49655   1.516  0.13089    
groupSM plot 2      195.54780    9.11190  21.461  < 2e-16 ***
groupSS plot 2      217.12857    9.11190  23.829  < 2e-16 ***
groupGE plot 3       21.68513    9.49655   2.283  0.02334 *  
groupSM plot 3      189.79613    9.49655  19.986  < 2e-16 ***
groupSS plot 3      192.53490    9.49655  20.274  < 2e-16 ***
time:groupGW plot 5  -0.01905    0.18816  -0.101  0.91944    
time:groupGW plot 6  -0.10231    0.18816  -0.544  0.58718    
time:groupGE plot 1  -0.80129    0.19357  -4.139 4.94e-05 ***
time:groupSM plot 1  -2.36256    0.18816 -12.556  < 2e-16 ***
time:groupSS plot 1  -2.55767    0.18816 -13.593  < 2e-16 ***
time:groupGE plot 2  -0.43459    0.19357  -2.245  0.02574 *  
time:groupSM plot 2  -2.34722    0.18816 -12.474  < 2e-16 ***
time:groupSS plot 2  -2.45552    0.18816 -13.050  < 2e-16 ***
time:groupGE plot 3  -0.56657    0.19357  -2.927  0.00378 ** 
time:groupSM plot 3  -1.82099    0.19357  -9.407  < 2e-16 ***
time:groupSS plot 3  -2.16847    0.19357 -11.202  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.51 on 223 degrees of freedom
Multiple R-squared:  0.9586,    Adjusted R-squared:  0.9544 
F-statistic: 224.7 on 23 and 223 DF,  p-value: < 2.2e-16
# Residual checking:
par(mfrow=c(2,2))
plot(M4) # a clear pattern in the residuals is observed (model assumptions are not met) 

# 6 - NULL MODEL (M5) ######################################################################################

M5 <- lm(lfmc ~ 1, data = lfmc)
summary(M5)

Call:
lm(formula = lfmc ~ 1, data = lfmc)

Residuals:
   Min     1Q Median     3Q    Max 
-88.75 -58.05 -31.45  52.66 231.06 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   95.645      4.919   19.45   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 77.3 on 246 degrees of freedom
IC_tab(M2, M3, M4, M5)#, weights = T, delta = TRUE, base=T, sort = TRUE) # the nonlinear mixed-effects model is clearly the best
IC_tab(M1)


par(mfrow=c(2,2), cex = 0.65) # (Figure 5)
v1 <- c(1, 2, 3, 4)
v2 <- c("GW", "GE", "SM", "SS")

# Nonlinear mixed-effects model (M1) 
rM1 <- resid(M1, type = "normalized")
fM1 <- fitted(M1)
plot(fM1, rM1, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)  
boxplot(rM1 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)
plot(rM1 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)
qqnorm(rM1, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM1)

# Nonlinear fixed-effects model (M2)
par(mfrow=c(2,2), cex = 0.65) # (Figure 5)

rM2 <- resid(M2)
fM2 <- fitted(M2)
plot(fM2, rM2, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)  
boxplot(rM2 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)
plot(rM2 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)
qqnorm(rM2, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM2)

par(mfrow=c(2,2), cex = 0.65)

# Linear mixed-effects model (M3) 
rM3 <- resid(M3, type = "normalized")
fM3 <- fitted(M3)
plot(fM3, rM3, xlab="Fitted LFMC (%)", ylab="Residuals", ylim=c(-3,3))
abline(a=0, b=0, lwd=2)  
boxplot(rM3 ~ lfmc$leaf.type, ylab="Residuals", xlab="Leaf type", xaxt="n", ylim=c(-3,3))
axis(side=1, at=v1, labels=v2, las=1)
plot(rM3 ~ lfmc$time, ylab="Residuals", xlab="Time (days)", ylim=c(-3,3))
abline(a=0,b=0, lw=2)
qqnorm(rM3, main="", ylab="Sample quantiles", xlab="Theoretical quantiles")
qqline(rM3)

par(mfrow=c(2,2), cex = 0.65)

# Clasical regression (M4) 
rM4 <- resid(M4)
fM4 <- fitted(M4)
plot(fM4, rM4, xlab="Fitted LFMC (%)", ylab="Residuals")
abline(a=0, b=0, lwd=2)  
boxplot(rM4 ~ lfmc$leaf.type, ylab="Residuals", xlab="", xaxt="n")
axis(side=1, at=v1, labels=v2, las=1)
plot(rM4 ~ lfmc$time, ylab="Residuals", xlab="Time (days)")
abline(a=0,b=0, lw=2)
qqnorm(rM4, main="", xlab="Theoretical quantiles", ylab="Sample quantiles")
qqline(rM4)

—-

LS0tCnRpdGxlOiAiTm9ubGluZWFyIFJlZ3Jlc3Npb24gKEFyY2hvbnRvdWxpcyBhbmQgTWlndWV6KSBwYXBlciIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKPGh0dHBzOi8vZmVtaWd1ZXouZ2l0aHViLmlvL25scmFhLWRvY3MvbmxyYWEtT2RkaS1MRk1DLmh0bWw+Cgo8aHR0cHM6Ly9vbmxpbmVsaWJyYXJ5LndpbGV5LmNvbS9kb2kvcGRmZGlyZWN0LzEwLjEwMDIvZWNlMy41NTQzPgoKCmBgYHtyfQojaW5zdGFsbC5wYWNrYWdlcygiZW1tZWFucyIpCgpsaWJyYXJ5KG5sbWUpCmxpYnJhcnkobmxyYWEpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShlbW1lYW5zKSAgICMgZW1tZWFuczo6Y29udHJhc3QoKSwgZW1tZWFuczo6ZW1tZWFucygpCmBgYAoKIyBTdGVwcyBpbiBmaXR0aW5nIGEgbm9ubGluZWFyIG1peGVkIG1vZGVsCgojIFBhcnQgSS4KCmBgYHtyfQpkYXRhKGxmbWMpCmxmbWMuZ2QgPC0gZ3JvdXBlZERhdGEobGZtYyB+IHRpbWUgfCBncm91cCwgZGF0YSA9IGxmbWMsIG9yZGVyLmdyb3VwcyA9IEZBTFNFKSAKbGZtYy5nZCRwbG90ID0gd2l0aChsZm1jLmdkLCBmYWN0b3IocGxvdCwgbGV2ZWxzID1jKCI0IiwgIjUiLCAiNiIsICIxIiwgIjIiLCAiMyIpKSkKaGVhZChsZm1jLmdkKQpgYGAKCmBgYHtyfQpnZ3Bsb3QobGZtYywgYWVzKHRpbWUsIGxmbWMpKSArCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG9lc3MiKSArCiAgZmFjZXRfd3JhcCh+bGVhZi50eXBlKSArCiAgeGxhYigiVGltZSAoZGF5cykiKSArCiAgeWxhYigiTE1GQyAoJSkiKQpgYGAKCiMjIE5PTkxJTkVBUiBNSVhFRC1FRkZFQ1RTIE1PREVMIChNMSkKCmBgYHtyfQpubHNMIDwtIG5sc0xpc3QobGZtYyB+IFNTZGxmKHRpbWUsIEEsIHcsIG0sIHMpLCBkYXRhID0gbGZtYy5nZCkKY29lZihubHNMKQpgYGAKCmBgYHtyLCBmaWcud2lkdGggPSA0LCBoZWlnaHQgPSAzfQpwbG90KGludGVydmFscyhubHNMKSkKYGBgCgojIyBSYW5kb20tRWZmZWN0cyBNb2RlCgpgYGB7cn0KbmwucmUgPC0gbmxtZShubHNMLCBjb250cm9sID0gbmxtZUNvbnRyb2wobWF4SXRlciA9IDUwMDAsIG1zTWF4SXRlciA9IDE1MDApKQpgYGAKCmBgYHtyfQpubC5yZS4xIDwtIHVwZGF0ZShubC5yZSwgcmFuZG9tID0gcGREaWFnKEEgKyB3ICsgbSArIHMgfiAxKSkKZnhmIDwtIGZpeGVmKG5sLnJlLjEpCmZ4ZgpgYGAKCmBgYHtyfQojIG5sLm1lIDwtIHVwZGF0ZShubC5yZS4xLCBmaXhlZCA9IGxpc3QoQSArIHcgKyBtICsgcyB+IGxlYWYudHlwZSksCiMgICAgICAgICAgICAgICAgIHN0YXJ0ID0gYyhBID0gZnhmWzFdLDAsMCwwLAojICAgICAgICAgICAgICAgICAgICAgICAgICAgdyA9IGZ4ZlsyXSwwLDAsMCwKIyAgICAgICAgICAgICAgICAgICAgICAgICAgIG0gPSBmeGZbM10sMCwwLDAsCiMgICAgICAgICAgICAgICAgICAgICAgICAgICBzID0gZnhmWzRdLDAsMCwwKSkKYGBgCgpgYGB7cn0KbmwucmUuMiA8LSB1cGRhdGUobmwucmUuMSwgcmFuZG9tID0gcGREaWFnKEEgKyB3ICsgbSB+IDEpKQpmeGYgPC0gZml4ZWYobmwucmUuMikKZnhmCmBgYAoKYGBge3J9Cm5sLm1lLjEgPC0gdXBkYXRlKAogIG5sLnJlLjIsIAogIGZpeGVkID0gbGlzdChBICsgdyArIG0gKyBzIH4gbGVhZi50eXBlKSwgIAogIHN0YXJ0ID0gYyhBPWZ4ZlsxXSwwLDAsMCwgdz1meGZbMl0sMCwwLDAsIG09ZnhmWzNdLDAsMCwwLCBzPWZ4Zls0XSwwLDAsMCkpCm5sLm1lLjEKYGBgCgpgYGB7cn0KbmwucmUuMgpgYGAKCmBgYHtyfQpBSUMobmwucmUuMiwgbmwubWUuMSkKYGBgCgpgYGB7cn0KSUNfdGFiKG5sLnJlLjIsIG5sLm1lLjEpCmBgYAoKYGBge3J9Cm5sLm1lLjIgPC0gdXBkYXRlKG5sLm1lLjEsIHJhbmRvbSA9IHBkRGlhZyhBICsgbSB+IDEpKQpJQ190YWIobmwubWUuMSwgbmwubWUuMikKYGBgCgpgYGB7cn0KbmwubWUuMyA8LSB1cGRhdGUobmwubWUuMiwgcmFuZG9tID0gcGREaWFnKEEgfiAxKSkgCklDX3RhYihubC5tZS4yLCBubC5tZS4zKQpgYGAKCmBgYHtyfQpwbG90KG5sLm1lLjMpCmBgYAoKYGBge3J9CnJzZSA8LSByZXNpZHVhbHMobmwubWUuMykgLyBzaWdtYShubC5tZS4zKQpoaXN0KHJlc2lkdWFscyhubC5tZS4zLCB0eXBlID0gIm5vcm1hbGl6ZWQiKSkKYGBgCgpgYGB7cn0KI3BhcihwdHkgPSAicyIpCnFxbm9ybShyc2UpCnFxbGluZShyc2UpCmBgYAoKYGBge3J9Cm5sLm1lLjMudm0gPC0gdXBkYXRlKG5sLm1lLjMsIHdlaWdodHMgPSB2YXJJZGVudChmb3JtID0gfjEgfCBsZWFmLnR5cGUpKQpJQ190YWIobmwubWUuMywgbmwubWUuMy52bSkKYGBgCgpgYGB7cn0KcGxvdChubC5tZS4zLnZtKQpgYGAKCmBgYHtyfQpoaXN0KHJlc2lkdWFscyhubC5tZS4zLnZtLCB0eXBlID0gIm5vcm1hbGl6ZWQiKSwgbWFpbiA9ICIiLCB4bGFiID0gIlJlc2lkdWFscyIpCmBgYAoKYGBge3J9CnFxbm9ybShyZXNpZChubC5tZS4zLnZtLCB0eXBlID0gIm5vcm1hbGl6ZWQiKSkKcXFsaW5lKHJlc2lkKG5sLm1lLjMudm0sIHR5cGUgPSAibm9ybWFsaXplZCIpKQpgYGAKCmBgYHtyfQpwbG90KEFDRihubC5tZS4zLnZtKSwgYWxwaGEgPSAwLjAxKQpgYGAKCiMjIEV2YWx1YXRpb24gb2YgRml4ZWQgRWZmZWN0cwoKYGBge3J9Cm5sLm1lLm1sIDwtIHVwZGF0ZShubC5tZS4zLnZtLCBtZXRob2QgPSAiTUwiKQpubC5tZS5tbApgYGAKCmBgYHtyfQpubC5tZS5tbC4xIDwtIHVwZGF0ZShubC5tZS5tbCwKICBmaXhlZCA9IGxpc3QoQSArIHcgKyBtIH4gbGVhZi50eXBlLCBzIH4gMSksCiAgc3RhcnQgPSBjKEEgPSBmeGZbMV0sIDAsIDAsIDAsIHcgPSBmeGZbMl0sIDAsIDAsIDAsIG0gPSBmeGZbM10sIDAsIDAsIDAsIHMgPSBmeGZbNF0pCikKbmwubWUubWwuMQpgYGAKCmBgYHtyfQpJQ190YWIobmwubWUubWwsIG5sLm1lLm1sLjEpIApgYGAKCmBgYHtyfQpubC5tZS5tbC4yIDwtIHVwZGF0ZShubC5tZS5tbCwgZml4ZWQgPSBsaXN0KEEgKyB3IH4gbGVhZi50eXBlLCBtICsgcyB+IDEpLCAgCiAgICAgICAgICAgICAgICAgICAgIHN0YXJ0ID0gYyhBPWZ4ZlsxXSwwLDAsMCwgdz1meGZbMl0sMCwwLDAsIG09ZnhmWzNdLCBzPWZ4Zls0XSkpIApubC5tZS5tbC4yCmBgYAoKYGBge3J9Cm5sLm1lLm1sLjMgPC0gdXBkYXRlKG5sLm1lLm1sLCBmaXhlZCA9IGxpc3QoQSB+IGxlYWYudHlwZSwgdyArIG0gKyBzICB+IDEpLCAgCiAgICAgICAgICAgICAgICAgICAgIHN0YXJ0ID0gYyhBPWZ4ZlsxXSwwLDAsMCwgdz1meGZbMl0sIG09ZnhmWzNdLCBzPWZ4Zls0XSkpIApubC5tZS5tbC4zCmBgYAoKYGBge3J9Cm5sLm1lLm1sLjQgPC0gdXBkYXRlKG5sLm1lLm1sLCBmaXhlZCA9IGxpc3QodyB+IGxlYWYudHlwZSwgQSArIG0gKyBzICB+IDEpLCAgCiAgICAgICAgICAgICAgICAgICAgIHN0YXJ0ID0gYyhBPWZ4ZlsxXSwgdz1meGZbMl0sMCwwLDAsIG09ZnhmWzNdLCBzPWZ4Zls0XSkpIApubC5tZS5tbC40CmBgYAoKYGBge3J9CklDX3RhYihubC5tZS5tbC4xLCBubC5tZS5tbC4yLCBubC5tZS5tbC4zLCBubC5tZS5tbC40KQpgYGAKCmBgYHtyfQpNMSA8LSB1cGRhdGUobmwubWUubWwuMiwgbWV0aG9kID0gIlJFTUwiKQpzdW1tYXJ5KE0xKQpgYGAKCmBgYHtyfQpmaXhlZihNMSkKYGBgCgpgYGB7cn0KcmFuZWYoTTEpCmBgYAoKYGBge3J9CmludGVydmFscyhNMSkKYGBgCgojIFBhcnQgSUkKCmBgYHtyfQojIDIuMi4xIC0gT3ZlcmFsbCBwcmVkaWN0aW9ucyAoVGFibGUgMykgLS0tLSAKQWd3IDwtIGZpeGVmKE0xKVsxXSAjICJBIiBmb3IgZ3Jhc3NlcyBpbiB0aGUgVyBzaXRlIChHVykgCkFnZSA8LSBmaXhlZihNMSlbMV0rZml4ZWYoTTEpWzJdICMgIkEiIGZvciBncmFzc2VzIGluIHRoZSBFIHNpdGUgKEdFKQpBc20gPC0gZml4ZWYoTTEpWzFdK2ZpeGVmKE0xKVszXSAjICJBIiBmb3IgTXVsbGludW0gc3Bpbm9zdW0gKFNNKQpBc3MgPC0gZml4ZWYoTTEpWzFdK2ZpeGVmKE0xKVs0XSAjICJBIiBmb3IgU2VuZWNpbyBmaWxhZ2lub2lkZXMgKFNTKQpXZ3cgPC0gZml4ZWYoTTEpWzVdICMgInciIGZvciBncmFzc2VzIGluIHRoZSBXIHNpdGUgKEdXKSAKV2dlIDwtIGZpeGVmKE0xKVs1XStmaXhlZihNMSlbNl0gIyAidyIgZm9yIGdyYXNzZXMgaW4gdGhlIEUgc2l0ZSAoR0UpIApXc20gPC0gZml4ZWYoTTEpWzVdK2ZpeGVmKE0xKVs3XSAjICJ3IiBmb3IgTXVsbGludW0gc3Bpbm9zdW0gKFNNKQpXc3MgPC0gZml4ZWYoTTEpWzVdK2ZpeGVmKE0xKVs4XSAjICJ3IiBmb3IgU2VuZWNpbyBmaWxhZ2lub2lkZXMgKFNTKQpNIDwtIGZpeGVmKE0xKVs5XSAjICJtIiBmb3IgYWxsIHRoZSBsZWFmIHR5cGVzClMgPC0gZml4ZWYoTTEpWzEwXSAjICJzIiBmb3IgYWxsIHRoZSBsZWFmIHR5cGVzCgpHVyA8LSBzdWJzZXQobGZtYywgbGVhZi50eXBlID09ICJHcmFzcyBXIikKR0UgPC0gc3Vic2V0KGxmbWMsIGxlYWYudHlwZSA9PSAiR3Jhc3MgRSIpClNNIDwtIHN1YnNldChsZm1jLCBsZWFmLnR5cGUgPT0gIk0uIHNwaW5vc3VtIikKU1MgPC0gc3Vic2V0KGxmbWMsIGxlYWYudHlwZSA9PSAiUy4gYnJhY3Rlb2xhY3R1cyIpCgpwYXIobWZyb3c9YygyLDIpLCBjZXg9MC43NSkKIyBHcmFzc2VzIFcKcGxvdChHVyR0aW1lLCBHVyRsZm1jLCB4bGFiPSIiLCB5bGFiPSJMRk1DICglKSIsIHlsaW09YygwLDEwMCksIG1haW49IkdyYXNzZXMgVyBzaXRlIiwgY29sPSJkYXJrZ3JlZW4iKQpjdXJ2ZSgoKEFndy1XZ3cpLygxICsgZXhwKChNLXgpL1MpKStXZ3cpLCBsd2Q9MiwgbHR5PTEsIGFkZD1ULCBjb2w9ImRhcmtncmVlbiIpIAojIEdyYXNzZXMgRQpwbG90KEdFJHRpbWUsIEdFJGxmbWMsIHhsYWI9IiIsIHlsYWI9IkxGTUMgKCUpIiwgeWxpbT1jKDAsMTAwKSwgbWFpbj0iR3Jhc3NlcyBFIHNpdGUiLCBjb2w9ImdyZWVuIikKY3VydmUoKChBZ2UtV2dlKS8oMSArIGV4cCgoTS14KS9TKSkrV2dlKSwgbHdkPTIsIGx0eT0xLCBhZGQ9VCwgY29sPSJncmVlbiIpCiMgTS4gU3Bpbm9zdW0KcGxvdChTTSR0aW1lLCBTTSRsZm1jLCB4bGFiPSJUaW1lIChkYXlzKSIsIHlsYWI9IkxGTUMgKCUpIiwgeWxpbT1jKDAsMzUwKSwgbWFpbj0iTS4gc3Bpbm9zdW0iLCBmb250Lm1haW49NCwgY29sPSJkYXJrb3JhbmdlIikKY3VydmUoKChBc20tV3NtKS8oMSArIGV4cCgoTS14KS9TKSkrV3NtKSwgbHdkPTIsIGx0eT0xLCBhZGQ9VCwgY29sPSJkYXJrb3JhbmdlIikgCiMgUy4gZmlsYWdpbm9pZGVzIApwbG90KFNTJHRpbWUsIFNTJGxmbWMsIHhsYWI9IlRpbWUgKGRheXMpIiwgeWxhYj0iTEZNQyAoJSkiLCB5bGltPWMoMCwzNTApLCBtYWluPSJTLiBicmFjdGVvbGFjdHVzIiwgZm9udC5tYWluPTQsIGNvbD0iZGFya3JlZCIpCmN1cnZlKCgoQXNzLVdzcykvKDEgKyBleHAoKE0teCkvUykpK1dzcyksIGx3ZD0yLCBsdHk9MSwgYWRkPVQsIGNvbD0iZGFya3JlZCIpIApgYGAKCmBgYHtyfQojIDIuMi4yIC0gUHJlZGljdGlvbnMgYXQgcGxvdCBsZXZlbCAoVGFibGUgMykgLS0tLQpBZ3cucDQgPC0gZml4ZWYoTTEpWzFdK3JhbmVmKE0xKVsxLDFdICMgIkEiIGZvciBncmFzc2VzIGluIHRoZSBwbG90IDQgKEdXKSAKQWd3LnA1IDwtIGZpeGVmKE0xKVsxXStyYW5lZihNMSlbMiwxXSAjICJBIiBmb3IgZ3Jhc3NlcyBpbiB0aGUgcGxvdCA1IChHVykgCkFndy5wNiA8LSBmaXhlZihNMSlbMV0rcmFuZWYoTTEpWzMsMV0gIyAiQSIgZm9yIGdyYXNzZXMgaW4gdGhlIHBsb3QgNiAoR1cpIApBZ2UucDEgPC0gZml4ZWYoTTEpWzFdK2ZpeGVmKE0xKVsyXStyYW5lZihNMSlbNCwxXSAjICJBIiBmb3IgZ3Jhc3NlcyBpbiB0aGUgcGxvdCAxIChHVykgCkFnZS5wMiA8LSBmaXhlZihNMSlbMV0rZml4ZWYoTTEpWzJdK3JhbmVmKE0xKVs3LDFdICMgIkEiIGZvciBncmFzc2VzIGluIHRoZSBwbG90IDIgKEdXKSAKQWdlLnAzIDwtIGZpeGVmKE0xKVsxXStmaXhlZihNMSlbMl0rcmFuZWYoTTEpWzEwLDFdICMgIkEiIGZvciBncmFzc2VzIGluIHRoZSBwbG90IDMgKEdXKSAKQXNtLnAxIDwtIGZpeGVmKE0xKVsxXStmaXhlZihNMSlbM10rcmFuZWYoTTEpWzUsMV0gIyAiQSIgZm9yIE0uIHNwaW5vc3VtIGluIHRoZSBwbG90IDEgKFNNKSAKQXNtLnAyIDwtIGZpeGVmKE0xKVsxXStmaXhlZihNMSlbM10rcmFuZWYoTTEpWzgsMV0gIyAiQSIgZm9yIE0uIHNwaW5vc3VtIGluIHRoZSBwbG90IDIgKFNNKQpBc20ucDMgPC0gZml4ZWYoTTEpWzFdK2ZpeGVmKE0xKVszXStyYW5lZihNMSlbMTEsMV0gIyAiQSIgZm9yIE0uIHNwaW5vc3VtIGluIHRoZSBwbG90IDMgKFNNKQpBc3MucDEgPC0gZml4ZWYoTTEpWzFdK2ZpeGVmKE0xKVs0XStyYW5lZihNMSlbNiwxXSAjICJBIiBmb3IgUy4gZmlsYWdpbm9pZGVzIGluIHRoZSBwbG90IDEgKFNNKSAgIApBc3MucDIgPC0gZml4ZWYoTTEpWzFdK2ZpeGVmKE0xKVs0XStyYW5lZihNMSlbOSwxXSAjICJBIiBmb3IgUy4gZmlsYWdpbm9pZGVzIGluIHRoZSBwbG90IDIgKFNNKSAKQXNzLnAzIDwtIGZpeGVmKE0xKVsxXStmaXhlZihNMSlbNF0rcmFuZWYoTTEpWzEyLDFdICMgIkEiIGZvciBTLiBmaWxhZ2lub2lkZXMgaW4gdGhlIHBsb3QgMyAoU00pCiAgCnBhcihtZnJvdz1jKDIsMiksIGNleD0wLjc1KSAjIChGaWd1cmUgNCkKIyBHcmFzc2VzIFcgc2l0ZSAKcGxvdChHVyR0aW1lLCBHVyRsZm1jLCB4bGFiPSIiLCB5bGFiPSJMRk1DICglKSIsIHlsaW09YygwLDEwMCksIG1haW49IkdyYXNzZXMgVyBzaXRlIiwgY29sPSJkYXJrZ3JlZW4iKQpjdXJ2ZSgoKEFndy5wNC1XZ3cpLygxICsgZXhwKChNLXgpL1MpKStXZ3cpLCBsd2Q9MSwgbHR5PTIsIGFkZD1ULCBjb2w9ImRhcmtncmVlbiIpICMgcGxvdCA0CmN1cnZlKCgoQWd3LnA1LVdndykvKDEgKyBleHAoKE0teCkvUykpK1dndyksIGx3ZD0xLCBsdHk9MiwgYWRkPVQsIGNvbD0iZGFya2dyZWVuIikgIyBwbG90IDUKY3VydmUoKChBZ3cucDYtV2d3KS8oMSArIGV4cCgoTS14KS9TKSkrV2d3KSwgbHdkPTEsIGx0eT0yLCBhZGQ9VCwgY29sPSJkYXJrZ3JlZW4iKSAjIHBsb3QgNgojIEdyYXNzZXMgRSBzaXRlCnBsb3QoR0UkdGltZSwgR0UkbGZtYywgeGxhYj0iIiwgeWxhYj0iTEZNQyAoJSkiLCB5bGltPWMoMCwxMDApLCBtYWluPSJHcmFzc2VzIEUgc2l0ZSIsIGNvbD0iZ3JlZW4iKQpjdXJ2ZSgoKEFnZS5wMS1XZ2UpLygxICsgZXhwKChNLXgpL1MpKStXZ2UpLCBsd2Q9MSwgbHR5PTIsIGFkZD1ULCBjb2w9ImdyZWVuIikgIyBwbG90IDEKY3VydmUoKChBZ2UucDItV2dlKS8oMSArIGV4cCgoTS14KS9TKSkrV2dlKSwgbHdkPTEsIGx0eT0yLCBhZGQ9VCwgY29sPSJncmVlbiIpICMgcGxvdCAyCmN1cnZlKCgoQWdlLnAzLVdnZSkvKDEgKyBleHAoKE0teCkvUykpK1dnZSksIGx3ZD0xLCBsdHk9MiwgYWRkPVQsIGNvbD0iZ3JlZW4iKSAjIHBsb3QgMwojIE0uIFNwaW5vc3VtIChFIHNpdGUpIApwbG90KFNNJHRpbWUsIFNNJGxmbWMsIHhsYWI9IlRpbWUgKGRheXMpIiwgeWxhYj0iTEZNQyAoJSkiLCB5bGltPWMoMCwzNTApLCBtYWluPSJNLiBzcGlub3N1bSIsIGZvbnQubWFpbj00LCBjb2w9ImRhcmtvcmFuZ2UiKQpjdXJ2ZSgoKEFzbS5wMS1Xc20pLygxICsgZXhwKChNLXgpL1MpKStXc20pLCBsd2Q9MSwgbHR5PTIsIGFkZD1ULCBjb2w9ImRhcmtvcmFuZ2UiKSAjIHBsb3QgMQpjdXJ2ZSgoKEFzbS5wMi1Xc20pLygxICsgZXhwKChNLXgpL1MpKStXc20pLCBsd2Q9MSwgbHR5PTIsIGFkZD1ULCBjb2w9ImRhcmtvcmFuZ2UiKSAjIHBsb3QgMgpjdXJ2ZSgoKEFzbS5wMy1Xc20pLygxICsgZXhwKChNLXgpL1MpKStXc20pLCBsd2Q9MSwgbHR5PTIsIGFkZD1ULCBjb2w9ImRhcmtvcmFuZ2UiKSAjIHBsb3QgMwojIFMuIGZpbGFnaW5vaWRlcyAoRSBzaXRlKQpwbG90KFNTJHRpbWUsIFNTJGxmbWMsIHhsYWI9IlRpbWUgKGRheXMpIiwgeWxhYj0iTEZNQyAoJSkiLCB5bGltPWMoMCwzNTApLCBtYWluPSJTLiBicmFjdGVvbGFjdHVzIiwgZm9udC5tYWluPTQsIGNvbD0iZGFya3JlZCIpCmN1cnZlKCgoQXNzLnAxLVdzcykvKDEgKyBleHAoKE0teCkvUykpK1dzcyksIGx3ZD0xLCBsdHk9MiwgYWRkPVQsIGNvbD0iZGFya3JlZCIpICMgcGxvdCAxCmN1cnZlKCgoQXNzLnAyLVdzcykvKDEgKyBleHAoKE0teCkvUykpK1dzcyksIGx3ZD0xLCBsdHk9MiwgYWRkPVQsIGNvbD0iZGFya3JlZCIpICMgcGxvdCAyCmN1cnZlKCgoQXNzLnAzLVdzcykvKDEgKyBleHAoKE0teCkvUykpK1dzcyksIGx3ZD0xLCBsdHk9MiwgYWRkPVQsIGNvbD0iZGFya3JlZCIpICMgcGxvdCAzCmBgYAoKYGBge3J9CiMgMi4yLjMgLSBEZXJpdmF0aXZlcyAoZHJ5aW5nIHNwZWVkKSAtLS0tIApsZm1jLkdXIDwtIGV4cHJlc3Npb24oKEFndy1XZ3cpLygxICsgZXhwKChNLXgpL1MpKStXZ3cpICMgTEZNQyBkeW5hbWljcyBmb3IgZ3Jhc3NlcyBpbiB0aGUgVyBzaXRlCmRzLkdXIDwtIGRlcml2KGxmbWMuR1csICJ4IiwgZnVuY3Rpb24uYXJnID0gVFJVRSkgIyBEcnlpbmcgc3BlZWQgZm9yIGdyYXNzZXMgaW4gdGhlIFcgc2l0ZQoKbGZtYy5HRSA8LSBleHByZXNzaW9uKChBZ2UtV2dlKS8oMSArIGV4cCgoTS14KS9TKSkrV2dlKSAjIExGTUMgZHluYW1pY3MgZm9yIGdyYXNzZXMgaW4gdGhlIEUgc2l0ZQpkcy5HRSA8LSBkZXJpdihsZm1jLkdFLCAieCIsIGZ1bmN0aW9uLmFyZyA9IFRSVUUpICMgRHJ5aW5nIHNwZWVkIGZvciBncmFzc2VzIGluIHRoZSBFIHNpdGUKCmxmbWMuU00gPC0gZXhwcmVzc2lvbigoQXNtLVdzbSkvKDEgKyBleHAoKE0teCkvUykpK1dzbSkgIyBMRk1DIGR5bmFtaWNzIGZvciBNLiBzcGlub3N1bQpkcy5TTSA8LSBkZXJpdihsZm1jLlNNLCAieCIsIGZ1bmN0aW9uLmFyZyA9IFRSVUUpICMgRHJ5aW5nIHNwZWVkIGZvciBNLiBzcGlub3N1bQoKbGZtYy5TUyA8LSBleHByZXNzaW9uKChBc3MtV3NzKS8oMSArIGV4cCgoTS14KS9TKSkrV3NzKSAjIExGTUMgZHluYW1pY3MgZm9yIFMuIGZpbGFnaW5vaWRlcwpkcy5TUyA8LSBkZXJpdihsZm1jLlNTLCAieCIsIGZ1bmN0aW9uLmFyZyA9IFRSVUUpICMgIyBEcnlpbmcgc3BlZWQgZm9yIFMuIGZpbGFnaW5vaWRlcwoKeHZlYyA8LSBzZXEoMCwgMTAwLCBsZW5ndGggPSAxMDAwKQp5LmRzLkdXIDwtIGRzLkdXKHh2ZWMpCnkuZHMuR0UgPC0gZHMuR0UoeHZlYykKeS5kcy5TTSA8LSBkcy5TTSh4dmVjKQp5LmRzLlNTIDwtIGRzLlNTKHh2ZWMpCgpwYXIobWZyb3c9YygyLDIpLCBjZXg9MC43NSkgIyAoRmlndXJlIDQpCnBsb3QoeHZlYywgeS5kcy5HVywgeGxpbT1jKDAsODApLCB5bGltPWMoMCw0KSwgCiAgICAgeWxhYj0iRHJ5aW5nIHNwZWVkIiwgeGxhYj0iVGltZSAoZGF5cykiLCB0eXBlPSJuIiwgbWFpbj0iR3Jhc3NlcyBpbiB0aGUgVyBzaXRlIikKbGluZXMoeHZlYywgLTEqKGF0dHIoeS5kcy5HVywgImdyYWQiKSksIGx3ZD0yLCBsdHk9MSwgY29sPSJkYXJrZ3JlZW4iKQpwbG90KHh2ZWMsIHkuZHMuR0UsIHhsaW09YygwLDgwKSwgeWxpbT1jKDAsNCksIAogICAgIHlsYWI9IkRyeWluZyBzcGVlZCIsIHhsYWI9IlRpbWUgKGRheXMpIiwgdHlwZT0ibiIsIG1haW49IkdyYXNzZXMgaW4gdGhlIEUgc2l0ZSIpCmxpbmVzKHh2ZWMsIC0xKihhdHRyKHkuZHMuR0UsICJncmFkIikpLCBsd2Q9MiwgbHR5PTEsIGNvbD0iZ3JlZW4iKQpwbG90KHh2ZWMsIHkuZHMuU00sIHhsaW09YygwLDgwKSwgeWxpbT1jKDAsNCksIAogICAgIHlsYWI9IkRyeWluZyBzcGVlZCIsIHhsYWI9IlRpbWUgKGRheXMpIiwgdHlwZT0ibiIsIG1haW49Ik0uIHNwaW5vc3VtIiwgZm9udC5tYWluPTQpCmxpbmVzKHh2ZWMsIC0xKihhdHRyKHkuZHMuU00sICJncmFkIikpLCBsd2Q9MiwgbHR5PTEsIGNvbD0iZGFya29yYW5nZSIpCnBsb3QoeHZlYywgeS5kcy5TUywgeGxpbT1jKDAsODApLCB5bGltPWMoMCw0KSwgCiAgICAgeWxhYj0iRHJ5aW5nIHNwZWVkIiwgeGxhYj0iVGltZSAoZGF5cykiLCB0eXBlPSJuIiwgbWFpbj0iUy4gZmlsYWdpbm9pZGVzIiwgZm9udC5tYWluPTQpCmxpbmVzKHh2ZWMsIC0xKihhdHRyKHkuZHMuU1MsICJncmFkIikpLCBsd2Q9MiwgbHR5PTEsIGNvbD0iZGFya3JlZCIpCmBgYAoKYGBge3J9CiMgMi4zIC0tLS0gVGVzdGluZyBkaWZmZXJlbmNlcyBhbW9uZyBsZWFmIHR5cGVzIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0gCgojIDIuMy4xIC0gUGFyYW1ldGVyICJBIiAtLS0tCmNvbnRyYXN0KGVtbWVhbnMoTTEsIH5sZWFmLnR5cGUsIHBhcmFtID0gIkEiKSwgInBhaXJ3aXNlIikKYGBgCgpgYGB7cn0KIyAyLjMuMiAtIFBhcmFtZXRlciAidyIgLS0tLSAKY29udHJhc3QoZW1tZWFucyhNMSwgfmxlYWYudHlwZSwgcGFyYW0gPSAidyIpLCAicGFpcndpc2UiKQpgYGAKCmBgYHtyfQoKIyAzLjEuMSAtIEJhc2Ugbm9ubGluZWFyIG1vZGVsIC0tLS0KbmwuZmUuMCA8LSBubHMobGZtYyB+ICgoQSAtIHcpIC8gKDEgKyBleHAoKG0gLSB0aW1lKS9zKSkpICsgdywgCiAgICAgICAgICAgICAgIHN0YXJ0ID0gYyhBPTE1LCBtPTMwLCBzPS0xNywgdz01KSwKICAgICAgICAgICAgICAgZGF0YSA9IGxmbWMpICMgaGVyZSwgdGltZSBpcyB0aGUgb25seSBwcmVkaWN0b3IgdmFyaWFibGUgCgpiIDwtIGNvZWYobmwuZmUuMCkgIyBjb2VmZmljaWVudHMgb2YgdGhlIGJhc2UgbW9kZWwKCiMgMy4xLjIgLSBMZWFmIHR5cGUgZml4ZWQgZWZmZWN0IC0tLS0KbmwuZmUuMSA8LSBubHMobGZtYyB+ICgoQVtsZWFmLnR5cGVdIC0gd1tsZWFmLnR5cGVdKS8oMSArIGV4cCgobSAtIHRpbWUpL3MpKSkgKyB3W2xlYWYudHlwZV0sIAogICAgICAgICAgICAgIHN0YXJ0ID0gbGlzdChBPXJlcChiWzFdLCA0KSwgbT0yNSwgcz0tMTAsIHc9cmVwKGJbNF0sIDQpKSwKICAgICAgICAgICAgICBkYXRhID0gbGZtYykgIyBhbmQgZml0IGlzIG1hZGUgdXNpbmcgdGhlIGJhc2UgbW9kZWwncyBjb2VmZmljaWVudHMgYXMgc3RhcnQgdmFsdWVzCgpiMSA8LSBjb2VmKG5sLmZlLjEpICMgY29lZmZpY2llbnRzIG9mIHRoZSBtb2RlbCB3aXRoIGxlYWYgdHlwZSBhbmQgdGltZSBhcyBwcmVkaWN0b3IgdmFyaWFibGVzCgojIDMuMS4zIC0gUGxvdCBmaXhlZCBlZmZlY3QgLS0tLSAKbmwuZmUuMiA8LSBubHMobGZtYyB+ICgoQVtncm91cF0gLSB3W2dyb3VwXSkvKDEgKyBleHAoKG0gLSB0aW1lKS9zKSkpICsgd1tncm91cF0sIAogICAgICAgICAgICAgICBzdGFydCA9IGxpc3QoQT1yZXAoYyhiMVsxXSxiMVsyXSxiMVszXSxiMVs0XSksMyksIG09MjUsIHM9LTEwLCB3PXJlcChjKGIxWzddLGIxWzhdLGIxWzldLGIxWzEwXSksMykpLAogICAgICAgICAgICAgICBkYXRhID0gbGZtYykgIyB0aGUgbW9kZWwgaXMgZml0IHVzaW5nICJiMSIgYXMgdGhlIHN0YXJ0IHZhbHVlcwoKSUNfdGFiKG5sLmZlLjAsIG5sLmZlLjEsIG5sLmZlLjIpCmBgYAoKYGBge3J9CiMgUmVzaWR1YWwgY2hlY2tpbmc6CnBhcihtZnJvdz1jKDEsMSkpCnBsb3QobmwuZmUuMikgIyBoZXRlcm9jZWRhc3RpY2l0eSBpbiB0aGUgcmVzaWR1YWxzCmBgYAoKYGBge3J9Cmhpc3QocmVzaWQobmwuZmUuMiksIG1haW49IiIsIHhsYWI9IlJlc2lkdWFscyIpICMgc2xpZ2h0bHkgc2tldyBkaXN0cmlidXRpb24gICAKYGBgCgpgYGB7cn0KcXFub3JtKHJlc2lkKG5sLmZlLjIpKQpxcWxpbmUocmVzaWQobmwuZmUuMikpCmBgYAoKYGBge3J9Ck0yIDwtIG5sLmZlLjIKc3VtbWFyeShNMikKYGBgCgpgYGB7cn0KIyA0LjEuMSAtIEJhc2UgbGluZWFyIG1peGVkLWVmZmVjdCBtb2RlbCAtLS0tCmwubWUgPC0gbG1lKGxmbWMgfiB0aW1lICogbGVhZi50eXBlLCByYW5kb20gPSB+MSB8IHBsb3QsIGRhdGEgPSBsZm1jKQoKIyBSZXNpZHVhbCBjaGVja2luZzoKcGxvdChsLm1lKSAjIGhldGVyb2NlZGFzdGljaXR5IGluIHRoZSByZXNpZHVhbHMKYGBgCgpgYGB7cn0KcGxvdChBQ0YobC5tZSwgcmVzVHlwZSA9ICJub3JtYWxpemVkIiwgbmEuYWN0aW9uPW5hLm9taXQpLCBhbHBoYT0wLjA1LCBncmlkPVRSVUUpICMgdGVtcG9yYWwgY29ycmVsYXRpb24gaXMgb2JzZXJ2ZWQgYXQgbGFnIDEKYGBgCgpgYGB7cn0KIyA0LjEuMiAtIFZhcmlhbmNlIG1vZGVsbGluZyAtLS0tIApsLm1lLnZtIDwtIHVwZGF0ZShsLm1lLCB3ZWlnaHRzID0gdmFySWRlbnQoZm9ybT1+IDEgfCBsZWFmLnR5cGUpKQogICAgCiMgUmVzaWR1YWwgY2hlY2tpbmc6CnBsb3QobC5tZS52bSkgIyB2YXJpYW5jZXMgdGhlcmUgc2VlbSB0byBiZSB3ZWxsIG1vZGVsZWQKYGBgCgpgYGB7cn0KQUlDKGwubWUsIGwubWUudm0pICMgbW9kZWxpbmcgdGhlIHZhcmlhbmNlIGltcHJvdmVzIHRoZSBtb2RlbCBmaXQgCmBgYAoKYGBge3J9CiMgNC4xLjMgLSBUZW1wb3JhbCBjb3JyZWxhdGlvbiBtb2RlbGxpbmcgLS0tLSAgCmwubWUudm0udG0gPC0gdXBkYXRlKGwubWUudm0sIGNvcnJlbGF0aW9uPWNvckFSTUEocD0xLHE9MSwgZm9ybT1+MSkpCnBsb3QoQUNGKGwubWUudm0udG0sIHJlc1R5cGUgPSAibm9ybWFsaXplZCIsIG5hLmFjdGlvbj1uYS5vbWl0KSwgYWxwaGE9MC4wNSwgZ3JpZD1UUlVFKSAjIHRoZSB0ZW1wb3JhbCBkZXBlbmRlbmNlIGlzIG1vZGVsZWQKYGBgCgpgYGB7cn0KQUlDKGwubWUudm0sIGwubWUudm0udG0pICMgbW9kZWxpbmcgdGhlIHRlbXBvcmFsIGNvcnJlbGF0aW9uIGltcHJvdmVzIHRoZSBtb2RlbCBmaXQKYGBgCgpgYGB7cn0KcGxvdChsLm1lLnZtLnRtKSAKYGBgCgpgYGB7cn0KcGxvdChBQ0YobC5tZS52bS50bSwgcmVzVHlwZSA9ICJub3JtYWxpemVkIiwgbmEuYWN0aW9uPW5hLm9taXQpLCBhbHBoYT0wLjA1LCBncmlkPVRSVUUpIApgYGAKCmBgYHtyfQpoaXN0KHJlc2lkKGwubWUudm0udG0sIHR5cGU9Im5vcm1hbGl6ZWQiKSwgbWFpbj0iIiwgeGxhYj0iUmVzaWR1YWxzIikgCmBgYAoKYGBge3J9CnFxbm9ybShyZXNpZChsLm1lLnZtLnRtLCB0eXBlPSJub3JtYWxpemVkIikpCnFxbGluZShyZXNpZChsLm1lLnZtLnRtLCB0eXBlPSJub3JtYWxpemVkIikpCmBgYAoKYGBge3J9CiMgNC4xLjQgLSBFdmFsdWF0aW9uIG9mIHRoZSBmaXhlZC1lZmZlY3RzIC0tLS0gCmwubWUubWwgPC0gdXBkYXRlKGwubWUudm0udG0sIG1ldGhvZCA9Ik1MIikgIyB0aGUgbW9kZWwgaXMgZml0dGVkIGJ5IE1heGltdW0gbGlrZWxpaG9vZCAKbC5tZS5tbC4xIDwtIHVwZGF0ZShsLm1lLm1sLCB+Li10aW1lOmxlYWYudHlwZSkgIyB0aGUgbW9kZWwgaXMgcmVkdWNlZCByZW1vdmluZyB0aGUgaW50ZXJhY3Rpb24gInRpbWUgeCBsZWFmIHR5cGUiIApJQ190YWIobC5tZS5tbCwgbC5tZS5tbC4xKSAjIHRoZSBBSUMgY29tcGFyaXNvbiBzdWdnZXN0cyB0aGUgaW50ZXJhY3Rpb24gdGVybSB0byBiZSBpbXBvcnRhbnQKYGBgCgpgYGB7cn0KTTMgPC0gbC5tZS52bS50bSAKc3VtbWFyeShNMykKYGBgCgpgYGB7cn0KTTQgPC0gbG0obGZtYyB+IHRpbWUgKiBncm91cCwgZGF0YSA9IGxmbWMpCnN1bW1hcnkoTTQpCmBgYAoKYGBge3J9CiMgUmVzaWR1YWwgY2hlY2tpbmc6CnBhcihtZnJvdz1jKDIsMikpCnBsb3QoTTQpICMgYSBjbGVhciBwYXR0ZXJuIGluIHRoZSByZXNpZHVhbHMgaXMgb2JzZXJ2ZWQgKG1vZGVsIGFzc3VtcHRpb25zIGFyZSBub3QgbWV0KSAKYGBgCgpgYGB7cn0KIyA2IC0gTlVMTCBNT0RFTCAoTTUpICMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjCgpNNSA8LSBsbShsZm1jIH4gMSwgZGF0YSA9IGxmbWMpCnN1bW1hcnkoTTUpCmBgYAoKYGBge3J9CklDX3RhYihNMiwgTTMsIE00LCBNNSkjLCB3ZWlnaHRzID0gVCwgZGVsdGEgPSBUUlVFLCBiYXNlPVQsIHNvcnQgPSBUUlVFKSAjIHRoZSBub25saW5lYXIgbWl4ZWQtZWZmZWN0cyBtb2RlbCBpcyBjbGVhcmx5IHRoZSBiZXN0CmBgYAoKYGBge3J9CklDX3RhYihNMSkKYGBgCgpgYGB7cn0KCgpwYXIobWZyb3c9YygyLDIpLCBjZXggPSAwLjY1KSAjIChGaWd1cmUgNSkKdjEgPC0gYygxLCAyLCAzLCA0KQp2MiA8LSBjKCJHVyIsICJHRSIsICJTTSIsICJTUyIpCgojIE5vbmxpbmVhciBtaXhlZC1lZmZlY3RzIG1vZGVsIChNMSkgCnJNMSA8LSByZXNpZChNMSwgdHlwZSA9ICJub3JtYWxpemVkIikKZk0xIDwtIGZpdHRlZChNMSkKcGxvdChmTTEsIHJNMSwgeGxhYj0iRml0dGVkIExGTUMgKCUpIiwgeWxhYj0iUmVzaWR1YWxzIiwgeWxpbT1jKC0zLDMpKQphYmxpbmUoYT0wLCBiPTAsIGx3ZD0yKSAgCmJveHBsb3Qock0xIH4gbGZtYyRsZWFmLnR5cGUsIHlsYWI9IlJlc2lkdWFscyIsIHhsYWI9IkxlYWYgdHlwZSIsIHhheHQ9Im4iLCB5bGltPWMoLTMsMykpCmF4aXMoc2lkZT0xLCBhdD12MSwgbGFiZWxzPXYyLCBsYXM9MSkKcGxvdChyTTEgfiBsZm1jJHRpbWUsIHlsYWI9IlJlc2lkdWFscyIsIHhsYWI9IlRpbWUgKGRheXMpIiwgeWxpbT1jKC0zLDMpKQphYmxpbmUoYT0wLGI9MCwgbHc9MikKcXFub3JtKHJNMSwgbWFpbj0iIiwgeGxhYj0iVGhlb3JldGljYWwgcXVhbnRpbGVzIiwgeWxhYj0iU2FtcGxlIHF1YW50aWxlcyIpCnFxbGluZShyTTEpCmBgYAoKYGBge3J9CiMgTm9ubGluZWFyIGZpeGVkLWVmZmVjdHMgbW9kZWwgKE0yKQpwYXIobWZyb3c9YygyLDIpLCBjZXggPSAwLjY1KSAjIChGaWd1cmUgNSkKCnJNMiA8LSByZXNpZChNMikKZk0yIDwtIGZpdHRlZChNMikKcGxvdChmTTIsIHJNMiwgeGxhYj0iRml0dGVkIExGTUMgKCUpIiwgeWxhYj0iUmVzaWR1YWxzIikKYWJsaW5lKGE9MCwgYj0wLCBsd2Q9MikgIApib3hwbG90KHJNMiB+IGxmbWMkbGVhZi50eXBlLCB5bGFiPSJSZXNpZHVhbHMiLCB4bGFiPSJMZWFmIHR5cGUiLCB4YXh0PSJuIikKYXhpcyhzaWRlPTEsIGF0PXYxLCBsYWJlbHM9djIsIGxhcz0xKQpwbG90KHJNMiB+IGxmbWMkdGltZSwgeWxhYj0iUmVzaWR1YWxzIiwgeGxhYj0iVGltZSAoZGF5cykiKQphYmxpbmUoYT0wLGI9MCwgbHc9MikKcXFub3JtKHJNMiwgbWFpbj0iIiwgeGxhYj0iVGhlb3JldGljYWwgcXVhbnRpbGVzIiwgeWxhYj0iU2FtcGxlIHF1YW50aWxlcyIpCnFxbGluZShyTTIpCmBgYAoKYGBge3J9CnBhcihtZnJvdz1jKDIsMiksIGNleCA9IDAuNjUpCgojIExpbmVhciBtaXhlZC1lZmZlY3RzIG1vZGVsIChNMykgCnJNMyA8LSByZXNpZChNMywgdHlwZSA9ICJub3JtYWxpemVkIikKZk0zIDwtIGZpdHRlZChNMykKcGxvdChmTTMsIHJNMywgeGxhYj0iRml0dGVkIExGTUMgKCUpIiwgeWxhYj0iUmVzaWR1YWxzIiwgeWxpbT1jKC0zLDMpKQphYmxpbmUoYT0wLCBiPTAsIGx3ZD0yKSAgCmJveHBsb3Qock0zIH4gbGZtYyRsZWFmLnR5cGUsIHlsYWI9IlJlc2lkdWFscyIsIHhsYWI9IkxlYWYgdHlwZSIsIHhheHQ9Im4iLCB5bGltPWMoLTMsMykpCmF4aXMoc2lkZT0xLCBhdD12MSwgbGFiZWxzPXYyLCBsYXM9MSkKcGxvdChyTTMgfiBsZm1jJHRpbWUsIHlsYWI9IlJlc2lkdWFscyIsIHhsYWI9IlRpbWUgKGRheXMpIiwgeWxpbT1jKC0zLDMpKQphYmxpbmUoYT0wLGI9MCwgbHc9MikKcXFub3JtKHJNMywgbWFpbj0iIiwgeWxhYj0iU2FtcGxlIHF1YW50aWxlcyIsIHhsYWI9IlRoZW9yZXRpY2FsIHF1YW50aWxlcyIpCnFxbGluZShyTTMpCmBgYAoKYGBge3J9CnBhcihtZnJvdz1jKDIsMiksIGNleCA9IDAuNjUpCgojIENsYXNpY2FsIHJlZ3Jlc3Npb24gKE00KSAKck00IDwtIHJlc2lkKE00KQpmTTQgPC0gZml0dGVkKE00KQpwbG90KGZNNCwgck00LCB4bGFiPSJGaXR0ZWQgTEZNQyAoJSkiLCB5bGFiPSJSZXNpZHVhbHMiKQphYmxpbmUoYT0wLCBiPTAsIGx3ZD0yKSAgCmJveHBsb3Qock00IH4gbGZtYyRsZWFmLnR5cGUsIHlsYWI9IlJlc2lkdWFscyIsIHhsYWI9IiIsIHhheHQ9Im4iKQpheGlzKHNpZGU9MSwgYXQ9djEsIGxhYmVscz12MiwgbGFzPTEpCnBsb3Qock00IH4gbGZtYyR0aW1lLCB5bGFiPSJSZXNpZHVhbHMiLCB4bGFiPSJUaW1lIChkYXlzKSIpCmFibGluZShhPTAsYj0wLCBsdz0yKQpxcW5vcm0ock00LCBtYWluPSIiLCB4bGFiPSJUaGVvcmV0aWNhbCBxdWFudGlsZXMiLCB5bGFiPSJTYW1wbGUgcXVhbnRpbGVzIikKcXFsaW5lKHJNNCkKYGBgCgojIC0tLS0K